## Dependent variable: Unemployment rate
## fm.unemp <- unemp ~ wells + cumucap
## jobloss.pool <- plm(fm.unemp,txPanel,model = "pooling")
## ## Fixed effects:
## ## One-way
## jobless.fe <- plm(fm.unemp,txPanel,model = "within")
## jobless.re <- plm(fm.unemp,txPanel,model = "random")
## ## Two-way
## jobless.twfe <- plm(fm.unemp,txPanel,model = "within", effect = "twoways")
## summary(jobless.twfe);
## jobless.twre <- plm(fm.unemp,txPanel,model = "random",random.method = "amemiya", effect = "twoways")
## summary(jobless.twre)

job.logtwts <- pgmm(formula = fm.gmmemp,data = txPanel, model = "twosteps", effect = "twoways")
summary(job.logtwts);

fm.gmmemp <- emp ~ lag(emp,1) + wells + cumucap|lag(log(emp),2:11)
## jobloss.logtwts <- pgmm(formula = fm.gmmunemp,data = txPanel, model = "twosteps", effect = "twoways")
## summary(jobloss.logtwts);

## Remove two outliers : Harris and Dallas
txDataR <- subset(txData, county != "Harris" & county != "Dallas")
txPanelR <- pdata.frame(txDataR, index = c("county", "year"), drop.index = TRUE, row.names = TRUE)

job.retwfe <- plm(fm.emp,txPanelR, model = "within", effect = "twoways")
summary(job.retwfe);

job.retwre <- plm(fm.emp,txPanelR,model = "random",random.method = "amemiya", effect = "twoways")
summary(job.retwre)

## For counties which have well drilling activities, analyse the impact of gas, oil and other wells.

txShaleWell <- subset(txData, wells != 0)
txShaleWell <- pdata.frame(txShaleWell, index = c("county", "year"), drop.index = TRUE, row.names = TRUE)
fm.emp.welltype <- emp ~ gas + oil

well.twfe <- plm(fm.emp.welltype,txShaleWell, model = "within", effect = "twoways")
summary(well.twfe);

well.re <- plm(fm.emp.welltype,txShaleWell, model = "random")
summary(well.re);


txShaleCty = txData
for (i in 1:254) {
  if (all(wellmat[i,] == 0)) {
   cty =  dimnames(wellmat) [[1]][i]
  txShaleCty = subset(txShaleCty, txShaleCty$county != cty)
 }
}
txShaleCty <- pdata.frame(txShaleCty, index = c("county", "year"), drop.index = TRUE, row.names = TRUE)

cty.twfe <- plm(fm.emp.welltype,txShaleCty, model = "within", effect = "twoways")
summary(cty.twfe);

cty.twre <- plm(fm.emp.welltype,txShaleCty, model = "random", effect = "twoways")
summary(cty.twre);

wt.twfe <- plm(fm.emp.welltype,txPanel, model = "within", effect = "twoways")
summary(wt.twfe);

wt.twre <- plm(fm.emp.welltype,txPanel, model = "random", effect = "twoways")
summary(wt.twre);

phtest(cty.twfe,cty.twre)
phtest(wt.twfe,wt.twre)


pay.logtwts <- pgmm(formula = fm.gmmwage,data = txPanel, model = "twosteps", effect = "twoways")
summary(pay.logtwts);
## Random effects:

## One-way effect, from Swamy and Arora (1972)
job.re <- plm(fm.emp,txPanel,model = "random")
summary(job.re)
## Two-way effects, from Amemiya (1971),
job.twre <- plm(fm.emp,txPanel,model = "random",random.method = "amemiya", effect = "twoways")
summary(job.twre)
attach(job.twre)
pdf("~/HOLA/research/green-jobs/figures/twre.pdf", height=6, width=6)
plot(txPanel$wells,txPanel$emp, main = "Random effects model")
abline(coefficients, lwd = 2, col = "blue")
mtext(bquote( emp == .(coefficients[2]) * txPanel$wells + .(coefficients[1])), side=3, line=0) 
dev.off()
detach()
pay.logtwfe <- plm(fm.logwage,txPanel,model = "within", effect = "twoways")
summary(pay.logtwfe);


## a = sapply(subset(txPair, shalecty == 0 & postshale == 0, select = emp), mean)
## b = sapply(subset(txPair, shalecty > 0 & postshale == 0, select = emp), mean)
## c = sapply(subset(txPair, shalecty == 0 & postshale > 0, select = emp), mean)
## d = sapply(subset(txPair, shalecty >0  & postshale > 0, select = emp), mean)

## d-c-(b-a)

## e = sapply(subset(txPair, post06 ==0 & shale == 0, select = lemp), mean)
## f = sapply(subset(txPair, post06 ==0 & shale == 1, select = log(emp)), mean)
## g = sapply(subset(txPair, post06 ==1 & shale == 0, select = log(emp)), mean)
## h = sapply(subset(txPair, post06 ==1 & shale == 1, select = log(emp)), mean)

## h-g-(f-e)
## txData$date = as.yearmon(txData$year)
## xyplot(wells~date|county,data = txData[txData$county == "Harris",])

wageR ~ -lag(wageR,1) - emp + wells - lag(wells,1) + cumucap.pre - lag(cumucap.pre,1) - shaleInd - lag(shaleInd,1) - newCap - lag(newCap,1) 
fm.logwage <- log(wageR) ~ - lag(log(wageR),1) - log(emp) + wells - lag(wells,1) + cumucap.pre - lag(cumucap.pre,1) - shaleInd - lag(shaleInd,1)
fm.gmmwage <- log(wageR) ~ + lag(log(wageR),1) - lag(log(emp)) + wells + lag(wells,1) + cumucap.pre + lag(cumucap.pre,1) - shaleInd - lag(shaleInd,1)|lag(log(wageR),2:11)

## pdf("~/HOLA/research/green-jobs/latex/figures/dataPair.pdf", height=6, width=8)
## pairs(txPanel)
## dev.off()


## Examine the dist. of the data
summary(txPanel);

## summary indicates the total variation of the variable and the share of this variation that is due to the individual and the time dimensions.
summary(txPanel$emp)
summary(txPanel$wells)
summary(txPanel$cumucap)
head(Within(txPanel$emp),15)
head(between(txPanel$emp),15)

## pbltest(fm.emp, txData, alternative = "onesided")


## For "short" panels with small T and large n
pwartest(fm.emp, txData)
pwfdtest(fm.emp, txData)
## Test for cross-sectional dependence
pcdtest(fm.emp, txData, model = "within", effect = "twoways")
## Robust cov matrix estimation
coeftest(job.twre,vcovHC)
pbltest(fm.wage, txData, alternative = "onesided")## For "short" panels with small T and large n
pwartest(fm.wage, txData)
pwfdtest(fm.wage, txData)
## Test for cross-sectional dependence
pcdtest(fm.wage, txData, model = "within", effect = "twoways")

## Robust cov matrix estimation
coeftest(pay.twre,vcovHC)

## Graphs of pair.shale counties

## xyplot(emp~year, data = txPair.Shale, groups = attr(txPair,"index")[,2], type = "o", auto.key=list(title="County", space = "right", cex=1.0))

## Results for the Kapoor et al. (2007)
errorre <- spml(formula = fm.emp, data = txData, index = "county",listw = wmatl,  model = "random", lag = F, spatial.error = "kkp")
summary(errorre)


## Spatial lag & error model, random effects
errlagre <- spml(formula = fm.emp, data = txData, index = "county", listw = wmatl,  model = "random", lag = TRUE, spatial.error = "kkp")
summary(errlagre)



## ## Spatial individual effects mu_i only, fixed effects
## spacefe <- spml(formula = fm.emp, data = txData, index = "county", listw = wmatl, model = "within", effect = "individual", method = "eigen", na.action = na.fail, quiet = TRUE, zero.policy = NULL, interval = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE)
## summary(spacefe)

## ## Time effects gamma_t only, fixed effects
## timefe <- spml(formula = fm.emp, data = txData, index = "county", listw = wmatl, model = "within", effect = "time", method = "eigen", na.action = na.fail, quiet = TRUE, zero.policy = NULL, interval = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE)
## summary(timefe)


## Two-way spatial lag & error model, fixed effects
errlagtwfe <- spml(formula = fm.wage, data = txData, index = "county", listw = wmatl, listw2 = wmatl, model = "within", lag = T, spatial.error = "b", effect = "twoways")

summary(errlagtwfe)


## ## Spatial error model, random effects



##Spatial error & lag  model, random effects
GM.full <- spgm(formula = fm.emp, data = txData, index = "county", listw = wmat, lag = TRUE, moments = "fullweights", model = "random", spatial.error = TRUE)

summary(GM.full)
## Spatial Hausman test#
testH <- sphtest(x = fm.emp, data = txData, index = "county",listw = wmatl, spatial.model = "error", method = "GM")
testH

mod1 <- spgm(formula = fm.emp, data = txData, index = "county", listw = wmat, lag = F, moments = "fullweights", model = "random", spatial.error = T)

mod2 <- spgm(formula = fm.emp, data = txData, index = "county", listw = wmat, lag = F, moments = "fullweights", model = "within", spatial.error = T)

test2 <- sphtest(x = mod1, x2 = mod2)
test2

## GMM Method
fm.gmmemp <- emp ~ lag(emp,1) + wells + cumuwells + newcap + cumucap|lag(emp,2:99)

job.logtwts <- pgmm(formula = fm.gmmemp,data = txPanel, model = "twosteps", effect = "twoways")

summary(job.logtwts);
